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ABSTRACT 

H2 formation is an important process in post-shock regions, since H2 is an 
r^ , active participant in the cooUng and shielding of the environment. The onset 

of H2 formation therefore has a strong effect on the temperature and chem- 



O ical evolution in the post shock regions. We recently developed a model for 

^ H2 formation on a graphite surface in warm conditions. The graphite surface 

H 

"^ acts as a model system for grains containing large areas of polycyclic aromatic 

I I hydrocarbon structures. Here this model is used to obtain a new description 

I of the H2 formation rate as a function of gas temperature that can be im- 

plemented in molecular shock models. The H2 formation rate is substantially 
\l higher at high gas temperatures as compared to the original implementation 

^-H of this rate in shock models, because of the introduction of H atoms which are 

chemically bonded to the grain (chemisorption). Since H2 plays such a key 



o 

^^ role in the cooling, the increased rate is found to have a substantial effect on 



the predicted line fluxes of an important coolant in dissociative shocks [O I] 
at 63.2 and 145.5 /xm. With the new model a better agreement between model 



J-j and observations is obtained. Since one of the goals of Herschel/PACS will be 

to observe these lines with higher spatial resolution and sensitivity than the 
former observations by ISO-LWS, this more accurate model is very timely to 
help with the interpretation of these future results. 
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1 INTRODUCTION 

H2 formation is an important process in post-shock regions, since H2 is an active participant 
in the coohng and shielding of the environment. The onset of H2 formation therefore has a 
strong effect on the temperature and chemical evolution in the post shock regions. This letter 
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aims to re-evaluate the estimations and assumptions concerning H2 reformation in molecular 
shock models. We mainly focus on jump-type (J-type) shocks. These can reach velocities 
of several hundred km s~^ causing an increase of the gas temperature > 10^ K. This will, 
in turn, result in the dissociation of molecular hydrogen, in contrast with continuous-type 
(C-type) shocks where H2 stays intact. The dissociation of H2 leads to the first cooling of the 
region ( Kwan||1977 London et al.||1977 ) to roughly 4300 K. Since all molecular material that 
is responsible for cooling under normal circumstances (and in C-type shocks), is dissociated, 
O cooling is the dominant cooling mechanism ( Flower et al.|2003 ). At some point, H2 starts to 
reform again on the surfaces of dust grains. Because of its high abundance, H2 cooling is more 
efficient than O cooling and it cools the gas phase until roughly 100 K when cooling by H2O 
takes over. O emission is classically used as a shock tracer, since it requires H2 dissociation 



and high temperatures, representative for J-type shocks (-Eup ~ 230 — 330 K; e.g., Nisini 



et al. 2002). Observations of O emission are typically complemented by H2 observations. 



which gives more lines with energy levels spread out over ~ 500-15 000 K. The combination 
of both observables can, in principle, provide a good handle on the physical conditions in 
shocks. However, accurate predictions and models are needed to translate these observations 
to physical conditions. 

Detailed molecular shock models can help us with this translation. They provide predic- 
tions of observables depending on a set of input parameters describing the physical condi- 
tions. Currently, most shock models include a number of "educated" estimations for parame- 
ter values. This work aims at minimising the uncertainty on at least one of these parameters: 
the reformation rate of H2 on grain surfaces. As explained above, the exact conditions under 
which H2 reforms is an important component in the evolutionary sequence and therefore 
in the interpretation of the observations. Often in gas phase models, the efficiency, i.e., the 
fraction of H atoms hitting a grain that leaves the surface as H2, is assumed to be unity 
and independent of temperature and pressure. Some further assumptions are made on the 
sticking of hydrogen atoms on the grains. 

Grains can be of silicate or carbonaceous nature and they can be covered by a mantle of 
ice. Here we focus on bare grains, since processing of ices in shocks, like inertial sputtering, 
has probably lead to the evaporation of the ice mantles. The chemical nature of the grain 
surface has an effect on the H2 formation rate in terms of both sticking and formation 
efficiency. We recently developed a model for molecular hydrogen formation on a graphite 
surface in warm conditions, such as shocks, PDRs and regions around AGNs. The graphite 
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surface acts as a model system for grains containing large areas of polycyclic aromatic 



hydrocarbon (PAH) structures (Pendleton & AUamandola 2002). This type of dust takes 



up a large fraction of the interstellar dust surface area and its interaction with hydrogen 
has been studied substantially in laboratory and theoretical studies. Our kinetic Monte 
Carlo model was tested against four independent experimental observations (H2 abstraction 
cross section, saturation coverage, surface structures, and desorption behaviour) and a good 
agreement was obtained ( Cuppen fc Hornekaer||2008 Gavardi et al.||2009 ). The present letter 
extends this model to shock conditions by determining the H2 formation efficiencies under 
these conditions. The improved rate is then introduced into a shock code to study its effect 
on the temperature and chemistry. The line fluxes of the [O I] at 63.2 and 145.5 fim are 
found to be strongly affected by the new H2 formation. In a comparative study between 



model and observations for HH54, Giannini et al. (2006) found an overestimation of the 



63.2 /im emission by at least a factor of 50. Here we will compare observations of HH46 to 
our new model and a much better agreement is obtained. 



2 THE MONTE CARLO METHOD 

For a detailed description of method and its parametrisation, we refer to|Cuppen &: Hornekaer 



(2008) and Gavardi et al. (2009). In short, the model uses a continuous-time, random-walk 



Monte Carlo scheme ( Chang et al.||2005 ) where hydrogen atoms can physisorb or chemisorb 
to the substrate which is represented by a honeycomb structure mimicking the graphite 
surface. The atoms can move across the surface and desorb from the surface. Physisorbed 
atoms can become chemisorbed after crossing a small barrier. The barriers and binding 



energies are obtained from theoretical methods. In accordance with experiments (Hornekaer 



et al. 2006), hydrogen atoms can cluster together on the surface. Some dimer structures. 



two hydrogen atoms chemisorbed to carbon atoms in close vicinity of each other, are found 
to be more stable than two, well-separated, monomers. The same holds for more complex 
structures. The model takes structures up to tetramers explicitly into account. 

Physisorbed atoms are found to play a key role in the formation of these complex struc- 
tures, even close to room temperature, where the lifetime of physisorbed atoms on the surface 



is short. Since their diffusion rate is high (Bonfanti et al. 2007), they scan a large part of 



the surface, find a chemisorbed H atom, and can form a dimer structure. 



Cuppen & Hornekaer (2008) describe a comparison between the model results and H2 
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Figure 1. The sticking probability on grapliite for physisorbed (circles and solid curve) and chemisorbed (dashed curves) H- 
atoms with different cheniisorption barriers. The sticking probability of H-atoms to water ice is plotted for comparison (Burke 
& HoUenbach 1983). 



abstraction experiments from Zecho et al. (2002). In order to obtain agreement, a mechanism 



was introduced where the excess energy of the chemisorbing H atom does not dissipate 
instantaneously, but can be used to cross, for instance, the H2 formation barrier. 

The Monte Carlo algorithm starts by the deposition of the first atom on the surface. 
The deposition rate is determined by the gas phase abundance of atomic hydrogen, 'n.(H), 
the sticking fraction to the grain, S, and the velocity of atoms in the gas phase 



^^H 



/o/eigas 



(1) 



t. 



with k the Boltzmann constant, Tgas the gas phase temperature of the incoming atom, and 
rriH the mass of an H atom. Then all possible events and the corresponding rates, Ri, of the 
deposited atom are considered depending on its environment and binding state. The time 
at which the atom will undergo an event is determined by 

-If (2) 

with tcurrcnt the currcut time, at which the first atom is deposited, and X a random number 
between and 1. This will be compared against the time of the second deposition. The first 
event will happen and new possible events are evaluated. This cycle is repeated until the 
end of the simulation. Event rates are determined by their barriers Ei 

E.; \ 



Ri = uexp 



T. 



(3) 



gram ^ 



with Tgrain the surface temperature of the grain and u the attempt frequency (lO^'^ Hz). 



An overview of all barriers that are included in the model is given in Cuppen fc Hornekser 



(2008). We have increased the physisorption desorption energy from 500 to 580 K to mimic 



the increased binding energies at step edges and other defects. 
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3 STICKING PROBABILITY 



The original implementation of the H2 formation rate in the model of Flower et al. ( 1985 ) 



assumes the efficiency to be unity independent of temperature and pressure. For the sticking 



probability, an expression describing the sticking of H atoms to water ice is used (Burke & 



Hollenbach||1983 ) . This expression is plotted as function of gas temperature in Fig. [T] (dash- 
dotted curve). However, grains are most likely not covered with water ice in post-shock 
regions. Inertial sputtering efficiently removes the ice mantles even for low velocity shocks 



(Guillet et al. 2009) 



Hydrogen atoms stick to graphite with a different efficiency at high gas phase tem- 
perature than to water. As explained earlier, hydrogen atoms can either be physisorbed 
or chemisorbed to the surface with a barrier between both potential wells, depending on 
the local environment. The trapping to both states has different probabilities and different 



dependencies on gas temperature. Sha et al. (2005) performed quantum studies on the trap 



ping of H atoms on graphite surface. They present the sticking probability after 1.0 ps as 
a function of incident energy. Here only the chemisorption channel was taken into account, 
and physisorption was ignored. To obtain a generalised expression that can be included in 
the Monte Carlo model, their results are described in terms of the product of a sticking 
probability and probability to cross the chemisorption barrier 

(4) 



S. 



exp 



1 + 5 X 10-2,/T,as + T^n + 1 X IQ-^^TL 



chein 

-gram i - -- -^ -gas 

with temperatures in K. Figure [l] plots this chemisorption sticking probability with the 
dashed curves for -Echcm = and -Echem = 0.15 eV. The first corresponds to sticking in a 
para-dimer configuration; the latter to sticking of monomers and corresponds to the work 



by Shaet ah (2005). 



The remaining non-sticking atoms still have an opportunity to become trapped in the 



physisorption well. To obtain this sticking probability we follow Burke & HoUenbach (1983), 



who used the so-called soft-cube method (Logan & Keck 1968). An incoming gas phase 
hydrogen atom is followed according to 

F 



dt^ mn 



(5) 



with Zgas the z-coordinate of the gas phase atom (normal to the surface) and F the force 
acting on this particle. The approached grain atom moves according to 
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Figure 2. The combined H-atom sticking and H2 formation efBciency as a function of gas phase temperature for four different 
grain temperatures. The contributions of physisorbed (top panel) and chemisorbed (middle panel) atoms to H2 formation are 
separated. The bottom panel gives the total of 5 X 77. The black dashed curve gives S X -q which was originally implemented in 
the shock model. 



d\. 



grain 



dt^ 



UJ Z. 



m 



gram 



(6) 



gram 



with u = kOu/h where 6]j is the Debye temperature of the sohd and h Planck's constant 
divided by 2tt. The gas phase atom starts from a large distance with a velocity following a 
Maxwell-Boltzmann distribution and the grain atom oscillates starting with a random phase 
and an amplitude of 



grain, max 



\ 



2kT, 



gram 



u^m 



(7) 



gram 



This method is found to give a reasonable agreement with experimental sticking probabilities 
of inert gas phase particles sticking to metal surfaces and we can therefore use it for the 
non-reactive trapping of hydrogen to graphite (physisorption) . To determine the force, F, 
we use an accurate H-coronene physisorption potential, which was fitted to ab initio data on 



MP2 level (Bonfanti et al. 2007) and which accurately reproduces experimental desorption 



data of H physisorbed to graphite. The sticking probability that is determined in this way 
is plotted by the open circles in Fig. [1} The solid line represents a fit to the data: 
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>S'phys ~ (1 + O-iyTgas + ^grain + Ct2^gas — C^S^gasj ' [°) 

with ai = 4.2 x lO^^ K^^/^ oa = 2.3 x IQ-^ K'^ and ag = 1.3 x 10"^ K'^. 



4 Ha FORMATION EFFICIENCY 

The combination of the sticking probabihty and the H2 formation efficiency, 77, is determined 
by Monte Carlo simulations. Simulation runs for a grid of different gas temperatures, H- 
atom gas abundances, and surface temperatures are performed. The shock model assumes 
the surface temperature to remain at a constant value of 15 K. The main heating mechanism 
for grains in shock regions is through UV photons. For the density and velocity regime that 
is covered here hardly any FUV photons are produced in the dissociation of the gas and the 



dust remains cold (HoUenbach & McKee 1989). Therefore simulations are performed at a 



constant low grain temperature, but with varying gas temperature, which affects the initial 
sticking to the surface. 

After a short initial equilibration time, a constant H2 formation efficiency is reached. This 
equilibration time scales with n(H). For n(H) = 10^ cm~^, it takes roughly 5 x 10^^ years to 
reach steady state. In comparison, n(H) in the post-shock gas is between 2 x 10^ — 2 x 10^ 
cm~'^ and the onset for H2 reformation occurs at longer timescales than the n(if )-dependent 
equilibration time (> 2x). Figure |2] shows the final S x i] as a function of gas temperature 
for different surface temperatures. The final values are obtained by averaging the steady 
state values obtained in five independent Monte Carlo simulation runs, with different initial 
seeds. Very little variation between the individual runs was found. Little to no dependence 
on n(H) was observed over the relevant range. 

As explained earlier, H2 can form from chemisorbed or physisorbed atoms. In Fig. |2| we 
distinguish between both mechanisms. The top panel shows the contribution to H2 formation 
by physisorbed atoms, the middle panel by chemisorbed atoms, and the bottom panel plots 
the sum of both contributions. The physisorbed mechanism is dominant at low gas and 
surface temperatures. At high gas temperatures most of hydrogen atoms scatter from the 
grain and for high grain temperatures the residence time of the atoms becomes shorter 
than the time between two depositions. Chemisorption requires the crossing of a barrier and 
is only efficient at higher temperatures. The dashed curves in Fig. |2] represents the original 
implementation of the product of the sticking probability and the efficiency (= 1), resembling 
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the physisorption contribution. Above Tgas = 500 K, the old and the new results start to 

deviate because of the introduction of the chemisorption mechanism in the new model. 

The solid curves in Fig. [2] give fits to the simulation data. These descriptions can be 

included in shock models to accurately describe the H2 formation rate. We find 

0.8 

for the physisorption contribution for Tgrain < 20 K and 
0.45 exp(-^) 

Q '^ ^ -'gas' ('\r\\ 

•Jchem^chcm — 1 i O v lO'^^T^ 

for the chemisorption contribution. These expressions are valid for 2 x 10'^ ^ n(H) ^ 2 x 10^ 
cm~^. 



5 SHOCK MODEL 

5.1 Model description and updates 



To test the impact of the above results, the molecular shock model described in Guillet 



et al. (2009) has been chosen as an application. This model has been chosen since it is 
a sophisticated ID model in which chemistry and cooling is treated in parallel with the 
magneto hydrodynamics. It is capable of treating both C-type and J-type shocks in a self- 
consistent manner. Grains are included as well as sputtering of the grain mantles. For the 
purposes of this paper only dissociative J-type shocks are considered, since H2 reformation 
is not at play in C-type shocks. The H2 reformation scheme has been adopted as follows: 

-Roid(H2) = ry/(H)n(grain)7rr|^^i„S'„ater (11) 

where /(H) the H-atom flux, n(grain) the number of grains, and 7rrgj.ain the geometric grain 
cross section. The efficiency is always assumed to be equal unity, and the sticking coefficient 
is given by the dash-dotted curve in Fig. [TJ The model has been updated to include the new 
results for H2 reformation 

Kcw(H2) = /(H)n(grain)7rr^^^i„(S'phys?7phys + 5'chem'7chem) (12) 

using Eqs. |9] and 10 



5.2 Results 

Figure [Sjplots the relative H2 abundance and the O cooling contribution for a grid of different 
shock simulations using different initial conditions in terms of velocity and density. The 
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Table 1. Predicted H2 1-0 ^{l) and [O I] 63.2 ^m and 145.5 fini line fluxes in the old and new version. 



riH V H2 1-0 5(1) [O I] 63.2 ^m [O I] 145.5 fim 

(cm-3) (km/s) (Wm-2sr-i) (Wm-^sr"!) (Wm-^sr"!) 

old new old new old new 

10^ 30 2.97(-8) 2.97(-8) 2.97(-9) 2.32(-9) 1.18(-10) 0.84(-10) 



40 1.66(-8) 2.40(-8) 6.76(-7) 4.42(-7) 2.38(-8 

50 1.43(-8) 1.72(-8) 6.47(-7) 4.23(-7) 2.08{-8 

10"' 30 8.72(-8) 1.55(-7) 3.00(-6) 1.76(-6) 8.44(-8 

40 1.07(-7) 1.40(-7) 1.57(-6) 0.62(-6) 4.12(-8 

50 1.13(-7) 9.48(-8) 1.18(-6) 0.45(-6) 2.96(-8 

10'"^ 30 3.74(-7) 5.88(-7) 7.16(-6) 4.03(-6) 1.94(-7 

40 4.06(-7) 5.67(-7) 3.41(-6) 0.95(-6) 9.07{-8 

50 4.02(-7) 2.47(-7) 2.08(-6) 0.42(-6) 5.41(-8 



1.46(-8) 
1.27(-8) 
4.89(-8) 
1.49(-8) 
1.00(-8) 
1.10(-7) 
2.49(-8) 
0.95(-8) 



a{—b) refers to a X 10 *> 

figure clearly shows that H2 is reformed at earlier times, especially for high densities and 
velocities. H2 takes therefore over as the main coolant earlier and at higher temperatures 
than previously seen. In general, this results in higher emission in high-f , high- J lines of H2 
at the cost of emission in low- J lines. The changes in f =1-0 S{1) emission, for instance, are 
small and sensitive to the initial physical conditions as is shown in Table [T] Changes within 
a factor of 1.8 are obtained, in both directions. The total H2 cooling summed over all lines 
is affected by less than 8%. Since H2 takes over as the main coolant earlier than in the older 
models, the timescale for, for example, oxygen cooling is shorter. This has the effect that 
the two [O I] lines at 63.2 and 145.5 /im have been over-estimated by factors of 1.5-5 in the 
older models for the same conditions as is indicated in Table [H 

As discussed above, we assume that the grain temperature does not change throughout 
the shock. Even if the dust would be heated to temperatures above 20 K, this will be 
at early times when H2 reformation occurs mainly through the chemisorption mechanism 
which is less temperature dependent. By the time that the physisorption mechanism starts 
to dominate, the dust has probably cooled down to lower temperatures. A grid of models 
with S'phys''7phys = shows indeed that the H2 production is only affected at late times. The 
[O I] fluxes are a factor of 0.33-0.87 lower, mainly because the cold post-shock gas is not 
fully converted to H2. 



Nisini et al. (2002) have measured the far-infrared spectra (45-197 jjm.) of 28 low- 
luminosity young embedded objects using the Infrared Space Observatory. Their sample 
includes the detection of ten 63.2 /im [O I] and six 145.5 /im [O I] lines. HH46 IRS is a 
typical source among these and we will use this for observational comparison. Assuming 
uniform emission, the observed emission of (5.25 ± 0.29) x 10^^^ Wm~^ of the 63.2 fim line 
corresponds to (3.49 ± 0.19) x 10~^ Wm~^sr~^ for the 80" aperture of the LWS receiver of 
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Figure 3. The relative abundance of H2 (left axis in black) and the O cooling (right axis in orange) for a grid of different 
velocities (30, 40, and 50 km s~^) and initial densities (10'', 10*, and 10^ cm^''). The old model results are represented by the 
dashed lines, the new by the solid lines. 



ISO, the emission of the 145.5 /im hne corresponds to (1.9 ± 0.5) x 10^^ Wm^^ sr^^. These 
values can be directly compared to the values in Table [ij With the old implementation of 
the H2 formation rate, almost all dissociative shocks {v > 25 kms^^) overproduce the [O I] 
63.2 /im and 145.5 yum emission with respect to the observed flux by at least factors of 18 
and 11, respectively {v = 50 kms^^ and n-n = 10^ cm^^). For higher densities discrepan- 
cies increase rapidly. For the new model, the predicted [O I] line fluxes are smaller and the 
density dependence is also reduced. The observed values are approached by factors of 12 
and 5, respectively, with very little dependence on the density. High velocities appear to 
be required to approach the observed values. In general, the agreement of the new models 
with the observations is better than for the old model, especially considering that the best 
agreement in the old model is only obtained for very speciflc conditions whereas in the new 



model a larger range of initial densities and velocities produces a good agreement. Garcia 



Lopez et al. (2010) reported the H2 1-0 S{1) line flux in six different knots of HH46. These 
observations are performed with a much smaller beam (6 times 0"3) than the ISO observa- 
tions. The total flux corresponds to (8.4 ± 0.1) x 10"'^ Wm~^ sr~^ which approached within 
a factor of 3.4 by the models for v = 50 kms~^ and wh = 10^ cm~^ (best flt for [O I]). As 
mentioned before, the H2 1-0 S{1) flux is rather sensitive to the physical conditions. For this 
particular combination of conditions the agreement with the model is less than with the old 
model. For lower initial velocity an improved agreement between model and observations is 
obtained. 
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6 CONCLUSIONS 

To summarise, we have introduced a new model for the H2 formation rate on carbonaceous 
surfaces as a function of gas temperature at a constant low grain temperature. This is ob- 
tained by a state-of-the-art model of H absorption on a graphite surface which includes 
both physisorption and chemisorption of complex structures. At low gas temperatures, phy- 
sisorbed atoms can stick to the surface whereas chemisorbed atoms do not have enough 
energy to overcome the barrier for sticking. The H2 formation rate is dominated by recombi- 
nation of physisorbed atoms in this regime. At high gas temperatures, most incoming atoms 
either scatter from the surface or become trapped in a chemisorption well. The H2 forma- 
tion rate is therefore mainly determined by chemisorbed atoms. This new model for the H2 
formation rate is used as input for a molecular shock model. The rate is substantially higher 
at high gas temperatures as compared to the original implementation of this rate in shock 
models. This increased rate has a substantial effect on the predicted line fluxes of [O I] at 
63.2 and 145.5 ;um and a better agreement between model and observations is obtained than 
with the previous model. One of the goals of Herschel/PACS will be to observe these lines 
with higher spatial resolution and sensitivity than the former observations by ISO-LWS. 
This more accurate model will help with the interpretation of these future results. 
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